# set working directory to source file location
# All models run using PanelMatch version 1.0.0
# uncomment if not installed in previous file
#devtools::install_version("PanelMatch", version = "1.0.0", repos = "http://cran.us.r-project.org")

#library(PanelMatch)


pacman::p_load(tidyverse, lfe, dplyr, ggplot2, ggpubr, stargazer, PanelMatch, purrr, cowplot, furrr)
#install.packages("sensitivitymw")
library(sensitivitymw)

# combined panel
load("../Data/combined_panel.RData")

combined_panel <- as.data.frame(combined_panel)
combined_panel$period <- as.integer(combined_panel$period)
combined_panel$log_rev <- log(combined_panel$revenue)

# check staff and bureau covariates are the same
cor(combined_panel$total_bur_filings, combined_panel$total_staff_filings)
combined_panel$log_contracts <- log(combined_panel$total_bur_filings)
combined_panel$log_lobby <- log(combined_panel$total_bur_lobbyists)

# create treatment indicator
combined_panel$treat <- ifelse(combined_panel$prev_exp == 1 | combined_panel$staff_exp == 1, 1, 0)
combined_panel$treat <- ifelse(is.na(combined_panel$treat)==T,0,combined_panel$treat)

# correlation between treatment and covariates
mod <- lm(treat ~ log_contracts, data = combined_panel)
summary(mod)

coef(mod)[2] * min(combined_panel$log_contracts)
coef(mod)[2] * max(combined_panel$log_contracts)

(coef(mod)[2] * quantile(combined_panel$log_contracts, prob = 0.75)) / (coef(mod)[2] * quantile(combined_panel$log_contracts, prob = 0.25))


coef(mod)[2] * (quantile(combined_panel$log_contracts, prob = 0.75)-
                quantile(combined_panel$log_contracts, prob = 0.25))

# create matched sets
comb <- PanelMatch(lag = 6, 
                   time.id = "period", unit.id = "reg_lobby_core_id", 
                   treatment = "treat", 
                   #refinement.method = "none", #alternative none
                   refinement.method = "mahalanobis", #alternative none
                   data = combined_panel, 
                   match.missing = T, 
                   covs.formula = ~ I(lag(log_contracts,1:4)) + I(lag(log_lobby,1:4)) + I(lag(log_rev,1:4)), 
                   qoi = "att" ,
                   outcome.var = "log_rev",
                   lead = 0:4, 
                   forbid.treatment.reversal = T)

overview <- print(comb$att)
overview[,2] <- as.numeric(as.character(overview[,2]))
overview$id <- 1:nrow(overview)

test <- comb$att

test2 <- as.numeric(as.character(test[[1]]))

test3 <- filter(combined_panel, period %in% (overview[1,2]-6):(overview[1,2]+4) &
                  reg_lobby_core_id %in% test2 | reg_lobby_core_id == overview[1,1])

weights <- data.frame(reg_lobby_core_id = as.numeric(as.character(comb$att[[1]])),
                      weights = attributes(comb$att[[1]])[[1]])
  

test3$treated <- NA
test3$treated <- ifelse(test3$reg_lobby_core_id == overview[1,1], 1, 0)
test3$et <- test3$period - overview[1,2] 
test3$post <- ifelse(test3$et == 0, 1, 0)

test3 <- left_join(test3, weights, by = "reg_lobby_core_id")

test3 <- test3 %>%
  group_by(reg_lobby_core_id) %>%
  arrange(reg_lobby_core_id, period) %>%
  mutate(lag_rev = dplyr::lag(log_rev, 1, order_by = period),
         change_rev = log_rev - lag_rev)


# compute a DiD by hand to check estimates
did <- test3 %>%
  ungroup() %>%
  filter(treated == 1 & et == 0) %>%
  select(change_rev) -   test3 %>%
  filter(treated == 0 & et == 0) %>%
  ungroup() %>%
  #summarise(av_change = mean(change_rev,na.rm=T)) %>%
  summarise(av_change = weighted.mean(x = change_rev,
                                      w = weights,
                                      na.rm=T)) %>%
  select(av_change) 





sens_df <- data.frame(treat = test3 %>%
             ungroup() %>%
             filter(treated == 1 & et == 0) %>%
             select(change_rev),
           control = test3 %>%
        filter(treated == 0 & et == 0) %>% ungroup() %>%
        select(change_rev) )

check <- senmwCI(sens_df, gamma = 1)

iter_sens <- function(sens_parm){
  this_est <- senmwCI(sens_df, gamma = sens_parm)$PointEstimate
  return(this_est)
}

base <- iter_sens(sens_parm = 1)[1]

sens <- map_df(seq(1, 3, by = 0.1), 
               ~iter_sens(.))

sens <- data.frame(sens, base, did)




####################
# put into a function to iterate
iter_did <- function(event_id){
  
  #crnt_id <- filter(overview, reg_lobby_core_id == id)
  crnt_id <- filter(overview, id == event_id)
  
  test <- comb$att
  
  test2 <- as.numeric(as.character(test[[event_id]]))
  
  test3 <- filter(combined_panel, period %in% (crnt_id[1,2]-6):(crnt_id[1,2]+4) &
                    reg_lobby_core_id %in% test2 | reg_lobby_core_id == crnt_id[1,1])
  
  weights <- data.frame(reg_lobby_core_id = as.numeric(as.character(comb$att[[event_id]])),
                        weights = attributes(comb$att[[event_id]])[[1]])
  
    test3$treated <- NA
  test3$treated <- ifelse(test3$reg_lobby_core_id == crnt_id[1,1], 1, 0)
  test3$et <- test3$period - crnt_id[1,2] 
  test3$post <- ifelse(test3$et == 0, 1, 0)
  
  test3 <- left_join(test3, weights, by = "reg_lobby_core_id")
  
    test3 <- test3 %>%
    group_by(reg_lobby_core_id) %>%
    arrange(reg_lobby_core_id, period) %>%
    mutate(lag_rev = dplyr::lag(log_rev, 1, order_by = period),
           change_rev = log_rev - lag_rev)
  
  did <- data.frame(treat = test3 %>%
    ungroup() %>%
    filter(treated == 1 & et == 0) %>%
    select(change_rev) ,   
    control = test3 %>%
    filter(treated == 0 & et == 0) %>%
    ungroup() %>%
    #summarise(av_change = mean(change_rev, na.rm=T)) %>%
    summarise(av_change = weighted.mean(x = change_rev,
                                          w = weights,
                                          na.rm=T)) %>%
    select(av_change) )
  
  return(did)

}    

#plan(multisession, workers = 6)
all_sens <- map_df(overview[,4],
               ~ iter_did(.))



  all_sens$did <- all_sens$change_rev - all_sens$av_change
  
  mean(all_sens$did, na.rm = T)
  
  run_sens <- function(sens_parm,df){
   this_est <- senmwCI(df, gamma = sens_parm,
                       method = "t")$PointEstimate
    return(this_est)
   }
  
  sens <- map_df(seq(1, 5, by = 0.1), 
                 ~run_sens(., df = all_sens$did))
  
  sens$gamma <- seq(1, 5, by = 0.1)
  0.4 - (1.645*0.12)
  0.19 - (1.645*0.12)
  0.34 - (1.645*0.112)
  
  
  ggplot(sens, aes(x = gamma, y = minimum)) +
    geom_line( lty = 3) +
    geom_line(aes(x = gamma, y = maximum),
              inherit.aes = F, lty = 3) +
    theme_classic() +
    geom_hline(yintercept = 0, lty = 2) +
    labs(x = "Gamma",
         y = "Bounded ATT") +
    scale_x_continuous(breaks = 1:5,
                       labels = as.character(c(
                         "1x\n(Equally Likely)",
                                      paste(2:4,"x",sep=""),
                                      "5x More\nLikely")))
  
  
    ggsave("../images/Figure5.pdf",
           width = 7, height =  5)
  
  
